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ABSTRACT 


A 

This report documents progress made in refining and improving the finite-volume explicit time 

-&&&> vJ vllt'i-'ilS.V.Q -Sr- \ (Vn^QCXiWO . 

marching method (1,2, and 3 ) during the time period from January to May 1986. The work is 
done under NASA grant NAG 3-593. Previously, extension had been made to the finite volume 
method to 

\/ improve the accuracy of the calculation of total pressure in inviscid flow^(l). 

2. extend the method to allow the calculation of laminar and turbulent boundary layers in internal 
flows^(2). vfrvJ'O 

3. improve the shock capturing properties of the method by introducing a Mach number de- 
pendent interpolation scheme for the pressure used in the calculating the density (3). 


The current work extends these developments by 

using the new pressure interpolation scheme in two dimensional viscous calculations^ 

2. including a more complete description of the viscous stresses^ 

3. introducing a criteria for the transverse upwind differencing which is a function of the ratio of 
transverse and streamwise mass fluxes^ 9 

4. allowing the calculation of internal flow where boundary layers are present on both wall of the 
duct. 


Specifically, this report is broken up into three sections. Section-l-discusses-in-detail the manner 
in which the viscous stresses are evaluated in the non-orthogonal, non-uniform grid. Section-2-inr 

- v ' ^ ■ife- C' ' &< ' rv 

vestigates the convergence and presents results for calculations of laminar flow in a converging duct. iwc 





Section 3 presents results for calculations of transonic turbulent flow in a converging-diverging 

fr 100 

nozzle; the results are compared with Sajben's measurements and calculations by other authors. 



SECTION 1 CALCULATION OF VISCOUS FORCES 


The momentum equation for unsteady flow in differential form is 

-^M.+ Fp -u U = - VP+ V-nVa+ V‘\i7u T + V • 5,,x4^- (1.1) 

ot J dXfc 

The normal stresses associated with second coefficient of viscosity, X, will be neglected in this 
analysis and in any subsequent calculations. We will be concerned in this scetion with how the 
viscous terms ( F • p Vu and F • jx Vu 7 ) are evaluated using the non-orthogonal, physical mesh 
system that is incorporated in the present method. 

In a control volume analysis of a flow field, we are interested in the actual forces which act upon 
the control surfaces and the components of these force in the coordinate directions ( x,y, and z ) 
rather than the derivatives of the shear stresses as seen in Eq. 1.1. To transform the governing 
equations from differential form (Eq. 1.1) into an integral form, we use Gauss' theorem. In terms 
of some arbitrary vector, g>, Gauss' theorem says, 

JJJ F • (g dVol = JJcj) • ndA (1.2) 

The viscous terms ( F *p Vu and F • p Vu 7 ) in Eq. 1.1, are converted from differential form to 
control volume form using Eq. 1 .2 and the result is, 

JJJ( F • p Vu + F • \iVu r )dVol = Jj(n Vu-dd + V~Vu r ’ d£) (1 .3) 

The current two-dimensional scheme uses control volumes which are made up of four straight line 
segments ( see Fig. 1.1 ). The surface integral in Eq. 1.3 is simplified into a summation over the 
four sides of the control volume. The integral over the surface in Eq. 1.3 can therefore be repres- 
ented as 
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Vu’dA + p Vu T • d£) = 


I {&K'V V U+ Ak'^Vu ) 

K=\ 


(1.4) 


Each area vector in Eq. 1.4 is assigned a magnitude equal to the area of the face in question and a 
direction which points in the outward normal direction. The evaluation of Eq. 1.4 results in the 
net viscous forces in the x and y directions for a control volume. 

For laminar flow, the absolute viscosity ,p, is used in Eq. 1.4 to evaluate the shear forces. For 
turbulent flow, an effective viscosity is used in Eq. 1.4. The effective viscosity that is used is the 
sum of the absolute viscosity and an apparent eddy viscosity. The current method uses a Prandtl 
mixing length model to evaluate this eddy viscosity. This mixing length model is outlined in Table 
1.1. The details of this mixing length calculation will be presented in section 1.3. 

The velocity gradients , Vu L , needed in Eq. 1.4 can be determined within a non-orthogonal grid by 
using, 


Vu, = — 


dUr 


Hk X Hi 


dui 


+ 


Hi x Hj 


dUr 


Hi ‘ (Hi x Hk) d/ Z2/ * {Hj x Hk) &J Z)/ • (Hj x Hk) 


(1.5) 


where D, t Dj , and D K are directional vectors along the grid directions (IJ, and K ), see Appendix 
A. For the two-dimensional case, duJdK = 0 and is a vector spanning the height of the duct 
( in the direction D, x Uj) Two typcial two dimensional control volumes are shown in Fig. 1.2. 
The directional vectors ( £), and Dj ) are identified in Fig. 1.2. The magnitudes of the vectors are 
dependent upon the grid spacing as can be seen in Fig. 1.2. The directional vectors that are shown 
in Fig. 1.2 would be used to calculate the velocity gradients applicable to the boundary common 
to both of the control volumes. 


The derivatives of the velocities in Eq. 1.5 are taken with respect to the grid indices, in other words, 
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Table 1.1 Prandlt Mixing Length Viscosity Model 


V-eff = + Ht 


u , » 


L is the smaller of 


0.08 times the width of the boundary layer 
0.41 times the distance to the nearest wall 
Van Driest Correction 

L = 0.41TO ~ expi - 126 ^]) 

Near Wall Correction 


( 1 ) 

( 2 ) 


( 3 ) 


V-eff = yjMVl + Hx) 


(4) 



(1.6) 


OUr 

—fij — ( u l)i+ 1 ~ ( u l)i 

where ( u L ), is the velocity at the beginning of vector £) , and (u L ) I+l is the velocity at the end of 
vector D,. For the geometries and boundary conditions investigated in the present work, the gra- 
dients in properties in the I direction are much smaller than gradients in the J direction. As a 
consequence, only the J derivative contributions to Eq. 1.5 will be considered in the present work. 

1.1 FORCES ON THE SOUTH FACE OF A CONTROL VOLUME 

In the actual calculations, for a non-uniform , non-orthogonal grid, the directional vectors are 
slightly different than those shown in Fig. 1.2. For the south face of a control volume, the area and 
directional vectors actually used are shown in Fig. 1.3. The directional vector ,Qj, and the velocity 
change are evaluated using the downstream nodal values because their use strengthens the 
centerpoint coefficient of the matrix of unknown variables. The directional vector J),, is located 
midway between the four nodes rather than at the boundary surface. This is because the viscosity 
and velocity gradient used in calculating the shear forces on the south face are evaluated midway 
between the two nodes. The shear stress is known to vary less through the boundary layer than the 
velocity gradient or the mixing length squared. Therefore, it is preferable to calculate the shear 
stress using a velocity gradient and mixing length midway between the grid points in the J-direction 
and then assign the resulting shear force to the face of the control volume between the points. The 
upper wall and lower wall control volumes are shown in Figs. 1.4 and 1.5 , respectively with the 
directional vectors identified. The shear stresses are evaluated midway between the wall and the 
near wall point and then the shear forces are assigned to the wall surface. The effective viscosity 
used to evaluate these wall shear forces is a combination of the laminar and turbulent viscosities 
given by 

V-eff= yjv-fol + Ft) (1-7) 
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This relationship is used only at the wall and has been shown ( Ref. 4 ) to allow a good calculation 
of wall shear stress with a near wall point further away from the wall than is typically required. 

Symmetry is used to calculate the forces on the north face of a control volume, in other words, 

*-► — ► 

FnorthJJ ~ ~ Fsouth,l,J+\ (1-8) 

where F narth and F south are the shear forces on the north and south faces of the control volume re- 
spectively. 

1.2 FORCES ON THE WEST FACES OF A CONTROL VOLUME 


The forces on the west face of the control volume are calculated by scalar multiplying an interpo- 
lated shear stress with the west face area vector ( see Eq. 1.4 ). Refering to Fig. 1.6, for a node point 
(I,J) which is located on the west face of control volume (I,J), first the velocity gradient and 
viscosity are calculated at the western edge of surfaces A and B which are midway between the node 
points (I,J) and (14-1) and node points (IJ + 1) and (14), respectively. The west sides of surfaces 
A and B are also the locations where the shear stresses are calculated for the north and south faces 
of the control volume (I- 1 ,J). The product of the velocity gradient and viscosity at node ( 14 ) is 
determined by linearly interpolating using the following interpolation formula, 


<Plj = <f>A + 


(y A /2) 


(i^ + ^L) 

' 2 2 ' 


(<P b ~ <Pa) 


(1.9) 


where <p A and <p B are the products , of the velocity gradient and the viscosity at the west sides of 
surfaces A and B respectively. Once these interpolated values have been calculated, Eq. 1.4 can 
be used to calculate the components of the shear forces on the west face. The shear stress is in- 
terpolated to node I,J rather than interpolating the velocity gradient since the shear stress varies less 
through the boundary layer than the velocity gradient. Similarly for the east face we have 
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( 1 . 10 ) 


FeasiJJ F\vest,I+ 1 J 

1.3 CALCULATION OF MIXING LENGTH AND VISCOSITY 

When the mixing length model of Table 1.1 is used in the 0.41y region, a distance normal to the 
wall, 'y"> must be determined to calculate the mixing length. For a flat wall, the distance to the 
wall from a point is, of course, measured along a line perpendicular to the wall. However, with the 
current grid system, grid lines and lines orthogonal to the wall are not coincident ( see Fig. 1.7 ). 
Fig. 1.7 shows 5 adjacent control volumes. The normal distance (L) between the north and south 
faces of a control volume is equal to the volume of the control volume divided by the distance (S) 
between the east and west faces of the control volume measured in the I grid direction. The total 
normal distance from the wall to the node point on the west face of a control volume ( ) is 

equal to the sum of all previous normal distances between that control volume and the nearest wall 
plus one half of the normal distance for that control volume. We can represent that distance as 

L we st= n i^k + 0.5 x (1.11) 

*=1 

where n is the number of control volumes from the wall. 

The total normal distance needed for calculating the mixing length for the south face of a control 
volume ( L south ) is equal to the sum of all previous normal distances to the previous node ( 
1 ) plus the equivalent normal distance half way between the two grid points in question. 
The procedure used to calculate this length is shown in Fig. 1.8. We can represent this total dis- 
tance as 


south/i 




est,n - 1 


(Ln - I + Ln) 


( 1 . 12 ) 
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To determine the mixing length in the outer part of the boundary layer, the boundary layer thick- 
ness measured normal to the wall must also be determined; The lengths that were calculated above 
can also be used in calculating the boundary layer thickness. The edge of the boundary layer is 
determined by using the magnitude of the normalized local total pressure gradient as a measure of 
its location. The normalized total pressure gradient used here is 


normalized total pressure gradient 


1 

"W 0.5p e u] 


(1.13) 


where "Ay" is the normal distance between two grid points. The local freestream values of density 
and velocity are used in Eq. 1.13. The grid line that is considered the freestream is input into the 
method. Starting from the freestream the normalized total pressure gradient ( Eq. 1.13 ) is calcu- 
lated between grid points and is compared with a characteristic normalized total pressure gradient. 
The edge of the boundary layer is defined as where the magnitude of the calculated gradient is equal 
to the characteristic gradient. The exact distance is determined using interpolation and is based 
upon the distances normal to the wall. The characteristic normalized total pressure gradient is de- 
termined by 


A^* t,char 1 _ 1 

*A y" 0.5p e Wg 

where CL is some characteristic length of the flow field, typically the duct width. 


(114) 


Once the mixing length has been calculated, the apparent eddy viscosity is calculated using Eq. 1 
in Table 1. The velocity gradient used in the eddy viscosity calculation is determined using the 
following expression, 


"du" _ ,A • Vu,2 , ,A- Vvj 

* ' J~UT ] ~1aT ] 


where A is the area vector on the south side of the control volume. This formula gives the mag- 
nitude of the total velocity gradient reflecting the non-orthogonality of the grid. 
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SECTION 2 LAMINAR BOUNDARY LAYER IN TWO 
CONVERGING DUCTS 


A la min ar boundary layer is calculated on the curved wall of two converging ducts. The other wall 
is straight and was treated as inviscid in the calculations. In all the calculations, the inlet boundary 
layer thickness is 5% of the inlet duct height. The inlet velocity profile is the Blasius profile. The 
inlet height of the duct is 44 mm and the exit height is 31 mm. The length of the duct is 180.4 
mm. The absolute viscosity is .001 kg/m s. The inlet freestream Mach number is 0.10. 

Two geometries are investigated with the basic dimensions given above. For one geometry, the 
curved wall radius is determined using a sine wave formulation. This geometry has a smooth 
transition from the inlet to the exit sections. The second geometry is essentially the same except 
that the radius is not determined using an analytical function and the second derivative of the wall 
radius is discontinuous. Both uniform and non-uniform grids are used in the boundary layer region. 
Figures 2.1 through 2.3 show the grid and geometry for the 3 arrangements to be investigated. The 
arrangements are 

1. Smooth geometry using non-uniform grid ( Fig. 2.1 ) 

2. Smooth geometry using uniform grid ( Fig. 2.2 ) 

3. Non-smooth geometry using uniform grid ( Fig. 2.3 ) 

A plot of the 2nd derivative of the wall radius as a function of axial position for the non-smooth 
and smooth geometries is seen in Fig. 2.4. The location of the discontinuity is also identified in 
Fig. 2.3. 

The specific ideas which are to be illustrated using this test case are, 
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1. Comparison of skin friction coefficients calculated using the finite volume method with those 
calculated using Thwaites method. 

2. Comparison of skin friction coefficients calculated using uniform and non-uniform grids. 

3. Comparison of skin friction coefficients calculated using smooth and non-smooth geometries. 

4. Investigation of pressure variations across the duct comparing the uniform and non-uniform 
grids results, and non-smooth geometry results. 

5. Investigation of inlet and exit boundary condition specifications. 

6. Investigation of the convergence properties of uniform vs non-uniform grids, coarse vs fine 
grids, and different multi-volume approaches. 

The radial distribution of grid points for the coarse, uniform and non-uniform grids and the fine 
non-uniform grid are tabulated in Table 2.1. 

For the smooth geometry in Fig. 2.1, the skin friction coefficient along the duct was calculated using 
the finite volume method and using the method of Thwaites (5). A coarse non-uniform grid was 
used for the finite volume calculations. Thwaites method is an integral method used to calculate 
the development of laminar boundary layers in incompressible flow. Using the wall static pressure 
distribution from the finite volume calculation as a boundary condition , the skin friction coefficient 
was calculated with Thwaites method. Fig. 2.5 compares the skin friction coefficients obtained 
from the finite volume calculations with those obtained using Thwiates method. The agreement 
is good for most of the flow ;however, near the exit of the duct there is a discrepancy. This dis- 
crepancy may be due to the fact that the boundary layer has grown to over 25% of the duct height 
at the exit. 
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Table 2.1 Radial Distribution of Grid Points 
for Laminar Calculations 



Uniform 

Grid 

Coarse 

Non-Uniform 

Grid 

Coarse 

Non-Uniform 

Grid 

Fine 


y/h 

y/h 

y/h 

PI 

.005 

.0024 

.0004 

P2 

.015 

.0071 

.00115 

P3 

.025 

.0142 

.00225 

P4 

.035 

.0283 

.0040 

P5 

.045 

.0563 

.0075 

P6 

.060 

.1125 

.0150 

P7 

.085 

.225 

.0300 

P8 

.125 

.400 

.0575 

P9 

.200 

.625 

.1125 

P10 

.375 

.875 

.225 

Pll 

.625 

- 

.400 

PI 2 

.875 

- 

.625 

P13 

— 

— 

.875 
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The skin friction coefficients calculated with the finite volume method using uniform and non- 
uniform grids are shown in Fig. 2.6. The agreement is excellent and it is concluded that grid spacing 
does not affect the prediction of the skin friction coefficient for this laminar flow. 


The skin friction coefficients for the flow through the smooth and non-smooth geometries are 
shown in Fig. 2.7. Superficially, the geometries in Figs. 2.2 and 2.3 look almost identical ; however, 
for the non-smooth geometry the second derivative of the lower wall radius of curvature is discon- 
tinuous (see Fig. 2.4) causing the solution to be erratic at the discontinuity. To ensure smooth 
calculations, geometries input into this program should have continuous second derivatives. 


The cause for the rapid increase in the skin friction coefficient at the discontinuity can be seen in 
Fig. 2.8. Fig. 2.8 is a plot of a transverse pressure coefficient which is defined as 


Cp t - 


P-Pfs 
0.5p u 2 


( 2 . 1 ) 


where p is the local static pressure at an axial location, p fs is the freestream static pressure at this 
axial location , and 0.5pw 2 is the dynamic pressure at the inlet to the nozzle. The different curves 
in Fig. 2.8 represent the pressure coefficient at every other grid line from the wall to the freestream. 
Because of curvature, the pressure gradient across the duct first rises and then falls. The wall 
pressure (pi ) is seen to change very rapidly at the point of the discontinuity. This large axial 
pressure gradient results in a large change in velocity at the discontinuity and thus the large change 
in skin friction coefficient that was seen in Fig. 2.7. 


For the smooth geometry the transverse pressure coefficient is shown in Fig. 2.9. The pressure 
changes are smooth throughout the entire length of the duct. The pressure coefficients shown in 
Fig. 2.9 are from calculations using a non-uniform grid. Fig. 2.10 shows identical results for cal- 
culations made using a uniform grid in the boundary layer. 
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In Figs, 2.8, 2.9, and 2.10, a small transverse pressure gradient can be observed at the inlet to the 
duct. This small transverse pressure gradient is a result of specifying the inlet v-velocity to be equal 
to zero. This boundary condition would be typical of inviscid calculations ; however, the boundary 
layer in this viscous flow does have a finite v-velocity as a result of the growth of the boundary 
layer. The growth of the boundary layer for this test case is large because of the high absolute 
viscosity used. An additional calculation was made by modifying the boundary condition on the 
inlet v-velocity. From a previous calculation, the flow angles at the second axial station were re- 
corded. These flow angles were then used as the inlet boundary condition for the inlet flow angle 
in a subsequent calculation. The transverse pressure coefficients for this case are shown in Fig. 2. 1 1. 
The inlet transverse pressure coefficient is reduced considerably when compared with the pressure 
coefficient shown in Fig. 2.9. When the Blasius solution was used to calculate the inlet v-velocity 
distribution, the solution over compensated and the transverse pressure gradient was in the opposite 
direction to that calculated when a zero v-velocity was specified. 

It was also found to be important that the inlet and exit lengths of the straight portions of the duct 
were long enough. Oscillations in the static pressure would occur at the inlet and exit if these 
lengths were too short. 

An additional pressure coefficient was calculated for each of the cases described above. The pres- 
sure coefficient is defined as 




= Po\~P 
0.5 pu 2 


( 2 . 2 ) 


where po i is the freestream total pressure, p is the local static pressure, and 0.5p u 2 is the inlet dy- 
namic pressure. This pressure coefficient again shows the effect of curvature on the pressure 
through the boundary layer but it also shows the relative acceleration of the flow through the 
nozzle. The plots are 


1. non- smooth geometry with uniform grid ( Fig. 2.12 ) 
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2. smooth geometry with non-uniform grid ( Fig. 2.13 ) 


3. smooth geometry with uniform grid ( Fig. 2.14 ) 

4. inlet flow angle specified with non-uniform grid ( Fig. 2.15 ) 

This test case also provided us with an opportunity to investigate the stability and convergence 
properties of various multi- volume approaches (2) within the boundary layer. The various strate- 
gies are 

1. [MV1] no multi- volume used in a uniform grid 

2. [MV2J multi- volume used in a uniform grid with time steps based upon properties for the in- 
dividual control volume farthest away from the wall. 

3. [MV3] multi-volume used in a uniform grid with time steps based upon average properties for 
the entire multi-volume region. 

4. [MV4J multi-volume used in a non-uniform grid with time steps based upon properties for the 
individual control volume farthest away from the wall. 

5. (MV 5] multi-volume used in a non-uniform grid with time steps based upon average properties 
for the entire multi- volume region. 

For the uniform and non-uniform grids, the multivolume approach was used for the 6 control 
volumes nearest to the wall. 

Two measures of convergence will be used in the following analysis. The maximum change 
in velocity in the flow field for one iteration divided by an average velocity for the flow field 
all multiplied by the time factor used in the time step determination will be our momentum 
residual, in other words, 
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X TIMEF 


(13) 


momentum residual = — — ^ max 

u avg 

The time factor, TIMEF, is included so that calculations using different time factors in the time 
step determination can be judged on a common basis. A time factor of 4.0 was used for all of 
the laminar calculations. The mass flow rate error is defined as 

mass flow error = | — ^-| max x 100 (2-4) 

The momentum residual for multi- volume approaches 1 through 4 are shown in Fig. 2.16. 
The mass flow error for approaches 1 through 4 are shown in Fig. 2.17. The non-uniform grid 
shows superior convergence properties when compared with the uniform grid results. When 
the uniform grid is used, approach MV3 is the best. 

When the momentum residual and mass flow error for the non-uniform grid approaches MV4 
and MV5 are compared in Figs. 2.18 and 2.19, the results are very similar. This agreement 
justifies the use of the simpler multi-volume approach ( MV4 ) for the non-uniform grid. The 
simpler multi-volume approach uses less computer time in the multi-volume part of the pro- 
gram. 

When a finer grid was used with the non-uniform grid arrangement in the boundary layer re- 
gion ( twice as many control volumes ), the skin friction coefficients were identical to those 
calculated from the coarse grid ( 5 points in the boundary layer). The momentum residual and 
the mass flow error for the coarse and fine grids are compared in Figs. 2.20 and 2.21 respec- 
tively. The convergence properties are very similar. 

The following conclusions can be reached from the above test case: 
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a. the geometries should be inspected for discontinuities in the second derivative to ensure 
smooth solutions. 

b. uniform and non-uniform grids give essentially identical results for this laminar boundary 
layer. 

c. the non-uniform grid (with a factor of 2 change in spacing) has superior convergence 
properties. 

d. the simple form of the multi- volume used with the non-uniform grid (MV4) is preferred 
because of the reduced computational effort required. 
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SECTION 3 SAJBEN S DIFFUSER CALCULATIONS 


Additional calculations have been made for Sajben's transonic diffuser [ 6,7 ]. Previous calculations 
using this geometry were presented in Refs. 1 and 2. These were preliminary calculations to show 
that the method could be used to calculate viscous transonic flow. Refinements and corrections 
have been made to the code since those progress reports and the following results reflect those re- 
visions and improvements. The improvements made to the method can be summarized as follows, 

1 . a new pressure interpolation scheme ( M&M formula ) is used to calculate the effective density 
( Ref. 3 ). 

2. a more complete treatment of the viscous stresses is included. 

3. the amount of upwinding used in the transverse direction ( transverse upwind differencing ) is 
based upon a new criteria which is a function of the ratio of mass fluxes through the transverse 
and streamwise faces of a control volume. 

4. an improved specification of the exit static pressure is used which reflects the side wall 
boundary layer blockage. 

5. the boundary layers are calculated on both walls of the diffuser. 

The diffuser geometry ( Model G ) is shown in Fig. 3.1; the throat height ,h, was 44 mm and the 
ratio of the exit height to throat height was 1.513. The calculations begin at x/h = -3.6 and end at 
x/h= 8.2. Fig. 3.1 also shows the computational grid used which had 87 grid points in the axial 
direction and 20 points across the flow. The radial distribution of grid points used is shown in 
Table 3.1. The development of a turbulent boundary layer was modeled on both the curved and 
the flat walls. The inlet boundary layer thicknesses were specified as 9 % and 4.5 % of the inlet 
diffuser height for the curved and flat wall boundary layers, respectively. An absolute viscosity of 
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Table 3.1 Radial Distribution of Grid Points 


1 

0.0005 

2 

0.0020 

3 

0.0050 

4 

0.0110 

5 

0.0230 

6 

0.0470 

7 

0.0950 

8 

0.1910 

9 

0.3080 

10 

0.4305 

11 

0.5695 

12 

0.6920 

13 

0.8090 

14 

0.9050 

15 

0.9530 

16 

0.9770 

17 

0.9890 

18 

0.9950 

19 

0.9980 

20 

0.9995 
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0.000018 kg/ m s was used. The inlet total pressure was 135 kPa and the inlet total temperature 
was 300 K. 


For this calculation, the ratio of the exit static pressure to the inlet total pressure was 0.826. This 
computational exit static pressure is equal to the experimental exit static pressure plus a correction 
in pressure for the side wall boundary layer blockage. Suction slots upstream of the throat and 
downstream of the exit plane reduce the side wall boundary layer blockage and improve the two- 
dimensionality of the flow. However, the side wall boundary layers do affect the effective flow area 
of the diffuser and the two dimensional computations must reflect this when a suitable exit static 
pressure is chosen. Therefore the exit static pressure of 0.826 is 1.5 % greater than the experimental 
static pressure because of the blockage effect of the side wall boundary layers. In the experiment, 
this test point results in transonic flow in the diverging portion of the duct with a Mach number 
of approximately 1.235 upstream of a nearly normal shock, and the flow remained fully-attached 
throughout the diffuser at this test condition. 

A plot of static pressure contours is shown in Fig. 3.2. The shock can be seen in the diverging 
portion of the duct. The shock is well defined as illustrated by the high clustering of contours at 
the shock. The M&M formula (3) is used to calculate the the interpolated pressure which is used 
in the calculation of the effective density. Fig. 3.3 shows a Mach number contour plot for the 
calculations. The extent of the boundary layer can be seen from this figure. 

The calculated and measured curved wall static pressures are compared in Fig. 3.4. The shock is 
very well defined and no overshoot occurs in the static pressure. The static pressures do not agree 
so well downstream of the shock because the 2-D calculations do not reflect the rapid increase in 
the side wall boundary layer thickness and its effective blockage. The point of minimum static 
pressure in the calculations is located at x/h= 1.5. This is taken to be the location of shock. The 
Mach number upstream of the shock was determined to be 1.256 from the calculated total pressure 
ratio across the shock in the freestream. 
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The computed and measured shock locations on the curved wall are compared in Fig. 3.5. The 
variable x au is the shock location on the curved wall and x gm is the shock location in the middle of 
the duct. The Mach number upstream of the shock, determined from static pressure measurements, 
is represented by M au in Fig. 3.5. 

The Mach number distribution through the nozzle at a fixed y/h of 0.0905 is shown in Fig. 3.6. 
This grid line was chosen because the maximum Mach number is located along it. The shock is 
sharp and no overshoots or undershoots occur. 

Comparisons of calculated and measured velocity profiles ( see Ref. 7) at four axial locations along 
the duct are shown in Figs. 3.7, 3.8, 3.9, and 3.10. The axial locations are x/h= 2.31, 4.03, 6.34, 
and 8.2 respectively. The agreement is good especially at the two downstream stations. The ve- 
locities in Figs. 3.7-3.10 are normalized with respect to the maximum computed or measured ve- 
locity in the duct at that axial location, which ever is applicable. In the calculations, the edge of 
the boundary layer is located where the normalized total pressure gradient is 25.0. It may be noted 
that the present calculations give much better agreeement with the measured velocity profiles than 
the calculations of Liou, Coakley, and Bergmann of Reference 8 ( see Fig. 3.17). Other investi- 
gators who have used Saj ben's diffuser as a test case have not shown a comparison of their com- 
puted velocity profiles with the measurements of Ref. 7. 

Fig. 3. 1 1 shows static pressure contours for calculations which were made using the three-point 
scheme for the interpolated pressure. The shock resolution is not nearly as well defined as that 
obtained when the M&M formula is used. Fig. 3.12 shows the corresponding Mach number con- 
tour plot when the three point interpolation scheme is used. 

The distribution of loss generation in the diffuser will be presented in three ways. First, the losses 
due to curved wall boundary layer, the flat wall boundary layer, and the shock loss in the freestream 
are compared in Fig. 3.13. These losses were determined by first calculating the mass flow rates in 
the boundary layers and the freestream at the diffuser exit using the calculated boundary layer 
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Table 3.2 


Example Calculation of Total Pressure 
Loss for a Boundary Layer 


HI £ 

5 (P . , -P , J • dm 
o t-mlet t-local 


m 


total 


* where is the total mass flow rate through 

a given cross-section of the diffuser, and m^ is the 
mass flow rate in the boundary layer at the exit. 
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thicknesses as the boundaries between regions. Then the mass averaged total pressure loss at each 
axial location was calculated by integrating the total pressure loss out to these fixed mass flow rate 
values and then normalizing these losses with respect to the inlet freestream total pressure ( see 
Table 3.2). All three losses are approximately the same. The curved wall ( bottom wall ) boundary 
layer contributes the largest proportion to the total losses because the boundary layer is thicker 
there. 

The total pressure loss along an inviscid streamline is compared with the mass averaged total pres- 
sure loss for the entire cross-section of the diffuser in Fig. 3.14. This figure allows another means 
of comparing the shock and total losses. The total pressure loss through the shock is very well 
defined. The mass averaged total pressure loss through the shock is approximately 30% greater 
than the shock loss alone. 

The mass averaged total pressure divided by the inlet freestream total pressure at the diffuser exit 
is calculated from the numerical results to be 0.9615. From data given to us by M. Sajben and T. 
Bogar, the mass averaged total pressure calculated from the experimentally measured data is 0.965. 
The experimental data used for this calculation was measured midway between the side walls. The 
agreement is good. 

The measured maximum Mach number in the exit plane of 0.51 agrees well with the calculated 
value of .511. The boundary layer thicknesses are measured to be approximately 25% of the duct 
height on the curved wall and 23 % of the duct height on the flat wall. The calculations determined 
the boundary layer thickness on the curved wall to be 25% of the duct height and 23% of the duct 
height of the flat wall. 

The above calculations are after 5000 iterations using a TIMEF of 4.0. It may have been possible 
to run the calculation with a TIMEF of 2.0 after the initial transients but because of computer cost 
this option was not attemped for this test case. The momentum residual for this calculation is 
presented in Fig. 3.15. The continuity error after 5000 iterations was less than .1%. The unusual 
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behavior observed between 3000-4000 iterations in Fig. 3. 1 5 appears to be typical of the method 
in general when transonic flow in calculated. The same general convergence behavior is seen for 
the 1-D nozzle ( Ref. 3 ) in Fig. 3.16. A time factor of 2.0 was used for these 1-D calculations but 
the same behavior is seen here on a smaller time scale. The total CPU time for the Sajben calcu- 
lations was approximately 35 minutes on the IBM 3031. 

A number of other workers have used this test case to verify the accuracy of their computational 
methods ( 8, 9, and 10 ). Liou, Coakley, and Bergmann's (8) calculations were made using a 
MacCormack type scheme in conjuction with a two-equation model for the turbulent stress pred- 
ictions. For the weak shock case, they used an exit static pressure of 0.8 x P tJnl „ . There was good 
agreement between their computed and measured wall static pressures. However, the calculated 
velocity profiles were not in good agreement with the experimental data. Fig. 3.17 shows a sample 
of their results. 

Liu, Shamroth, and McDonald (9) used this test case to verify their computational method which 
solves the governing equations using a consistently split linearized block implicit scheme. A mixing 
length model was used to model turbulent stresses. For the weak shock case, an exit static pressure 
of 0.807 x P tJnl „ was used. They compared their computed wall static pressures with the exper- 
imental pressures and found good agreement (see Fig. 3.18). The variable ,a, referred to in Fig. 3.18 
is the value of the artificial dissipation parameter used in their calculations. They found that for 
values of a less than 0. 1 the solution was invariant. Also presented in Ref. 9 were static pressure 
and Mach number contour plots for this weak shock case. These contour plots are similar to those 
in Figs. 3.2 and 3.3. The shock definition in the present work is sharper than Liu et. al. but fewer 
axial grid points were used in their calculations. No comparison was made by Liu, Shamroth, and 
McDonald between the measured and computed velocity profiles. 

Talcott and Kumar (10) also used the weak shock case for Sajben's diffuser as a test of their com- 
puter program's accuracy. They used an implicit MacCormack scheme to solve the governing 
equations and a Baldwin- Lomax mixing length turbulence model to predict the turbulent stresses. 
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They also used an exit static pressure of 0.807 * Pw Fig. 3.19 shows a sample of their results. 
They present a static pressure contour plot for this weak shock case. As with Liu et. al., the shock 
definition of Talcott and Kumar's method is not as sharp as that obtained using the present method 
but again they have used fewer axial grid points in their calculations. Also, the static pressure 
contours at the throat of the diffuser are erratic. A plot of wall static pressure is also presented. 
No comparision was made between the computed and measured velocity profiles. 

In summary, the following observations can be made about the present calculations of the weak 
shock case in Sajben's diffuser. 

1. There is better agreement between measured and calculated velocity profiles when the present 
method and static pressure specification are used. 

2. The present method compares favorably with other methods in the areas of static pressure and 
Mach number contour, predictions and wall static pressure prediction. 

3. The agreement between calculated and measured mass averaged total pressure losses for the 
diffuser is excellent. 
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Fig. 3.1 Geometry and Grid for Sajben's Diffuser Calculations 
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Fig. 3.17 SamDle Results from Liou, Coakley, and Bergmann (8) 
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Fig. 3.19 Sample Results from Talcott and Kumar (10) 
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APPENDIX A 


To calculate the velocity gradient, Vu L , for a nonorthogonal grid system, we must be able to 
transform gradients in the grid directions (I,J, and K) into gradients in the coordinate directions 
(x,y, and z). From the transformation law of a vector, the gradient of a scalar in one coordinate 
system can be related to the gradient in another coordinate system using the direction cosines be- 
tween the two coordinate systems. We will use the cartesian system as our base coordinate system 
here. Let Hijh, and D K be directional vectors which have a direction parallel with the non- 
orthogonal coordinate system and span a unit change in I,J, or K respectively. These vectors are 
shown in Fig. Al. They can be described in terms of their cartesian components. In other words, 
we can represent them as 


Dj = a x i + b^l + c x k 

(A. 1) 

Dj = a 2 i + l> 2 l+ c 2 k 

(A.2) 

Hk = a 3 i + 63 / + c 3 fc 

(A3) 


where i,J, and k are unit vectors in the cartesian coordinate directions, and a,b, and c are the 
magnitudes of the components of the directional vectors. 


The change of a variable ,u, from one point to another within a grid system can be represented as 

■fjT-K/+i -u r = jl +] dr- Vu (A. 4) 

where r is a position vector between the two nodes. Rather than dealing with a local derivative in 
a continuum , we are dealing with approximations of derivatives using finite changes. If the direc- 
tional vector, Qi, is chosen such that its termini correspond to the nodes at which the properties 
are evaluated in Eq. A.4, and taking Vu as uniform between grid points, 

A. 1 



= (£/+) - C/)* Vu. 


M.5) 


But Hi = C/+i - & therefore, 


# =a- r “ 

04-6) 


04-7) 

ir - Fu 

(A.8) 


where I,J, and K are the coordinate indices for the non-orthogonal coordinate system. Expanding 
equations A6,A7, and A8, we get 


5m _ „ 5m , » 5m , „ 5m 

TT ~ 37 + *'17 + c '-37 


5m _ „ 5m , » 5m , , 5m 

a7 " 37 + ^ 02 37 


5m _ „ 5m , > 5m , „ 5m 

IF " + + C3 1F 


04-9) 
(A 10) 
0411) 


where a, = £),^, the x component of D,\ b t = D,^ etc... Equations A9,A10, and All can be com- 
bined into a single matrix equation, 


[ du' 


\ b x q 

[ du 

51 



5xr j 

5u 


a 2 h c 2 

‘ 

/5m 

\ 5J 




5j; 

du 


a 3 &3 c 3 


5m 

dK) 


- - 


5z ' 


A. 2 
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For the calculation of viscous stresses, Eq. 1.5, we need the gradients of velocity ,Vu L , in cartesian 
coordinates. To obtain this gradient, we need to take the inverse of the above matrix [A], in other 
words, 


fM.) 


b x 

c \ 


( 8u \ 

dx 






81 

8u 

l = 

a 2 

b2 

C 2 


8u 

8y 

\ 





8J 

8u 


* 3 

63 

C3 

1 

8u 

8z) 


— 


A 

l 

8K 
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The inverse of matrix A , [A]~ l , will now be determined. From linear algebra, 


[A]’ 1 =-L-(adj\A)) 


(A. 14) 


where \A 1 is the determinant of [A] and adj[/l] is the adjoint of [ A]. The determinant of [/l] is 


detM] = 


a \ b \ c, 

<h h c 2 

a 3 b 3 c 3 


*1 


bi c 2 

h c 3 


“ b \ 


c 2 
<*3 c 3 


+ C X 


a 2 b 2 

*3 b 3 


a l (^2 C 3 ~ hCl) ~ b l( a 2 c 3 ~ «3«2) + C 1 ( a 2 b 3 ~ a 3h) 


(A. 15) 


This can be identified as 


det[i4] — D i ' ( Dj x Dr) 


(AA6) 


Now let us represent the adjoint of A , adj [A], in the form 


adj[A] = 


A u A\2 A\2 

A 21 A 22 A 23 
A31 A 32 A 33 


(AM) 
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where the components of the adjoint matrix are 


A\\ “ 


*2 c 2 
h c 3 


= ^c 3 - b 3 c 2 


(A. 18) 


However, An, can be identified as the x-component of 22/ x R K , in other words, 


62 C 3 - b 3 c 2 = x component of Rj x Rk 


(A. 19) 


The components of the adjoint matrix, /4 12 and /4 !3 ,are the y and z components of Rj x Z)*. Both 
the results from the numerator and denominator of Eq. A. 12, lead to the result that 


Rj * Rk d u L + Rk x Ri d u L + Ri * Rj d u L 
Ri • ( Rj x Rk) Rj • (Rj x 22 a') SJ 22 / * (22/ x 22^) dK 


(A. 20) 


This result can be generalized to the other two components of velocity, v and w. This is the result 
that is used to calculate velocity gradients in the current program. 
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